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ABSTRACT 

We study the heating of the cool cores in galaxy clusters by cosmic-rays (CRs) acceler- 
ated by the central active galactic nuclei (AGNs). We especially focus on the stability 
of the heating. The CRs stream with Alfven waves in the intracluster medium (ICM) 
and heat the ICM. First, assuming that the heating and radiative cooling is balanced, 
we search steady state solutions for the ICM and CR profiles of clusters by solving a 
boundary value problem. The boundary conditions arc set so that the solutions are 
consistent with observations of clusters. We find steady state solutions if the magnetic 
fields are strong enough and the association between the magnetic fields and the ICM 
is relatively weak. Then, we analyse the stability of the solutions via a Lagrangian 
perturbation analysis and find that the solutions are globally stable. We confirm the 
results by numerical simulations. Using the steady state solutions as the initial con- 
ditions, we follow the evolution of the profiles for 100 Gyr. We find that the profiles 
do not evolve on time scales much larger than cluster lifetimes. These results, as well 
as consistency with observations of radio mini-halos, suggest that the CR heating is a 
promising mechanism to solve the so-called "cooling flow problem" . 

Key words: cosmic rays — cooling flows — galaxies: clusters: general — galaxies: 
clusters: individual: A1795, A2052, A2199, A2597 



1 INTRODUCTION 

The radiative cooling time of the intracluster medium (ICM) 
in the cores of gal axy clusters is often smaller than the 
age of the clusters (Sarazi ni l 19861 ). If there were no heat- 
ing sources, massive cooling flo ws toward th e cluster cen- 
tres would develop in the cores l|Fabianll 19941 ). However, X- 
ray observations have revealed that such flows do not exist 
in the cores (e .g. Ikebe et al. 19971: Makishima et al. 200 ll: 



Cosmic-rays (CRs) accelerated by the AGNs 
are another promising carrie r of t he energy to the 
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Rephaeli fc Silkl 



Tucke r & Rosncr 1983; Rephaej Il987l: 



1995 



Colafrancesco. Par, fc Pe Ruiulal 



20041 ; IPfrommer et all 120071 ; lJubelgas et al l T2OO8) 



200ll ; ^Matsushi ta et al 



Peterson etHkl l200ll; iTamura et all l200ll ; iKaastra et~ai1 
2002). This means that some un- 



streaming i s one way to 
the ICM dRephaelil Il979t 



deposit their energ y 
Boehringer fc Morfil 



CR 
into 



1988; 



known heating sources prevent such flows from developing 
(cooling flow problem). 

AGNs at the cluster centres are often thought to be the 
heating sources. However, the energy transfer mechanism 
from the AGNs to the surrounding ICM is not understood. 
Conventional mechanical heating ( e.g. shocks, or sound 
waves) may cause thermal instability (iFuiita fc Suzukil2005l ; 
iMathews. Faltenbacher. fc Brighentil 20061 ) . Thus, strong 
turbulence would be required to stabilise the heating, if the 
cores are mainly mechanically heated. 
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Loewenstein. Zweibel fc Begelmanl 1 199 it iGuo fc Ohlbood ; 
Pfrommerll2013l ). Each CR particle moves with a velocity 
close to the light velocity c. However, in a cluster, they 
are scattered by Alfven waves in the ICM. Thus, the 
CRs effectively move along with the waves with a bulk 
(streaming) velocity v st , which is much smaller than c. Since 
the CR pressure is high around the central AGN and waves 
are excited in the direction of the CR pressure gradient (e.g. 
lLongairl Il994l ). the waves and the CRs scattered by them 
move outward in the cluster in general. As a result, the CR 
pressure does PdV work against the ICM, which ultimately 
heats the ICM. Using numerical simulations, we showed 
that the ICM can be stably heated by the CR streaming 
(|Fuiita fc Ohirall201ll . hereafter Paper I). The mam reason 
of the stability is that the CRs stream in the ICM and 
the heating is not localised around the AGN. In successive 
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studies, we calculated non-thermal emissions from the CRs 
l|Fuiita fc Ohiral |2012| . hereafter Paper II), and we showed 
that the radial profiles of radio mini-hal os observed in 
clust ers can be reproduced in our model l|Fuiita fc Ohiral 
l2013l . hereafter Paper III). 

In this paper, we study the reason of the stability of the 
CR heating in more detail based on a perturbation analy- 
sis and numerical simulations. This paper is organised as 
follows. In § [21 we explain our models and show the ba- 
sic equations. In § [3] we solve those equations and obtain 
steady state solutions when the CR streaming velocity is 
the Alfven velocity. Moreover, we analyse the stability of the 
steady state solutions. In § 21 we discuss the evolution of the 
ICM with time and radio mini-halos observed around cluster 
cores. §[5] is devoted to conclusions. Throughout this paper 
we assume a ACDM cosmology with Q, m — 0.3, Q,a = 0.7, 
and h = 0.7, where Ho = 100 h km s _1 Mpc -1 . We consider 
protons as CRs unless otherwise mentioned. 



2 MODELS 
2.1 Equations 

The model we adopt in this study is basically the same as 
that in Papers I— III. Assuming that the cluster is spherically 
symmetric, the flow equations are 



reason, we assume a simple cooling function: 
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where p is the gas density, u is the gas velocity, P g is the gas 
pressure, P c is the CR pressure, Pb is the magnetic pressure, 
g is the gravitational acceleration, k(T) is the coefficient for 
thermal conduction, T is the temperature, n e is the electron 
density, A is the cooling function, H a t is the heating by CR 
streaming, -ff co n is the heating by Coulomb and hadronic col- 
lisions, u is the CR transport velocity, D(p) is the diffusion 
coefficient for CRs averaged over the CR spectrum, Pi oss is 
the energy loss by Coulomb and hadronic collisions, and S c 
is the source term of CRs. Energy densities of the gas and 
the CRs are respectively denned as e 3 = Pg/ifg — 1) and 
e c = P c /(7c — 1), where 7 S (= 5/3) and 7 C are the adiabatic 
indices for the ICM and CRs, respectively. For the latter, we 
assume 7 C = 5/3, because we found that the energy spec- 
trum of the CRs must be steep and most of the CRs have 
low energies by comparing model predictions with radio ob- 
servations (Papers II and III). 

Since we later perform a perturbation analysis, we need 
to make the model simpler. We do not include thermal con- 
duction and CR diffusion (k — and D — 0). For the same 



n*A(T) = pC 



21 (If < 5 > 



ilRvbicki fc Lightmanl Il979l : iKim fc Naravanl 120031 ; 
iGuo. Oh. fc Ruszkowskill2008l ). We also ignore the gradient 
of magnetic pressure (9Pb /dr = 0) in equation ^ , because 
it is not dynamically important. 



2.2 AGN 

In our model, CRs are accelerated by the central AGN. The 
source term of the CRs is given by 

3 — V Z/AGN 



5* c = 



v -r\ \r J 



47T rf(ri/r ) 
x(l-e- (r/ro)2 )e- (r/ri) 



(6) 



where Lagn is the energy injection rate from the AGN. We 
assume that 



Lagn = —eMc 



2.5 x 10- 4 / v ioM 0y r 



-M 



7 '(7) 



where e is the parameter, and M is the flow rate of the gas. 
In this model, accelerated CRs are first carried by buoyant 
bubb les from the AGN out to large distances l|Guo fc Ohl 
2008). As the bubbles filled with the CRs adiabatically rise, 
the CRs may escape from the bubbles into the ICM or they 
may be injected into the ICM through the shredding of the 
bubbles by Rayleigh- Taylor and Kelvin-Helmholtz instabil- 
ities. Thus, the CRs are gradually injected into the ICM as 
the buoyant bubbles rise. Unless otherwise mentioned, we 
fixed the parameters at the values similar to the ones used 
in Papers I— III (e = 2.5 x 10" 4 , v = 3.2, r = 20 kpc, and 
n = 150 kpc), because they give results that are consis- 
tent with observations (see later). In Papers I and II, we 
assumed that the CRs are accelerated at the forward shock 
of a cocoon created through the AGN activities. However, 
it may be more appropriate to assume that CR protons are 
accelerated around the central black hole and that AGN 
jets or winds consist of the CR proto ns, especially when 
e is large (jj 3.1 of Pap er III; see also ISikora et all 120051 : 



e is large (s 6.L 01 rap 
iToma fc Takaharal f2012) 



2.3 Cosmic-rays 

The CR transport velocity is given by « = u + v st in 
equation (j4|. The most simple idea is that v st is the 
Alfven velocity va, because the CRs are scattered by the 
Alfven waves and move with them. However, there have 
been debates about this issue. It has been indicated that 
v st may be much larger than va, because in hot ICM, 
Alfven waves may damp at small wave lengths via inter- 
actions with thermal protons. In this case, the sound ve- 
locity of the gas c 3 may be more plau s ible a s the stream- 
ing v elocity l|Holman. Ionson. fc Scottl 1 19791 ; lEnfilin et al.l 
1201 if ) . However, there also have been indications that the 
deficit of the scattering at short wave lengths is over- 
come by other effects such as mirror scattering or wave 
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cascading (e.g. Felice fc Kulsrudl l200ll ; ISchlickeiserl |2002; 
I Wiener. Oh. fc Gud 20131 ) . If this is the case, the assump- 
tion of v st — va is appropriate. We consider the latter case 
(v st = va) in the following sections. We discuss the former 
case (v Bt = c s ) in Appendix. 

In this study, we investigate steady state solutions. 
Therefore, we assume that the growth of the waves is bal- 
anced with non-linear wave damping. In other words, the 
CR energy consumed to grow the waves is equal to the wave 
energy that is put into the ICM through the wave damping. 
Thus, the heating term of the CR streaming is given by 

(8) 



H at = — Vat- 



dr 



Note that v st > and dP c /dr < in our calculations. 

In Papers I and II, we studied the evolution of wave 
amplification. In those studies, we conservatively assumed 
that the balance between the wave growth and the damping 
is archived when the wave energy density reaches to that 
of the background magnetic fields, because some of non- 
linear damping mechanisms should be effective after that. 
As shown in Paper I (its Fig. 3), the time scale in which the 
growth and the damping is balanced is ~ Gyr. Since it is 
less than the time scale of a cluster lifetime (~ 10 Gyr), the 
assumption of the balance in this study is justified. More- 
over, if the damping is very efficient, the saturation may 
be achieved even faster. For example, if the time scale of 
the damping is proportional to the inverse of the gyro fre- 
quency of CR particle s (e.g. non-linear Landau damping; 
iFelice fc Kulsrudll200ll ). the saturated wave energy density 
is much smaller than that of the background magnetic fields. 

We ignore Coulomb and hadronic collisions (H co \\=0 
and Tioss = 0) in equations Q and (J4j , because they do not 
much affect the results as follows (see also Papers I and III). 
The collisional heating term is given by H co \\ — — F c — 1\/6, 
and the loss term is given by ri oss = — T c — Fh, where F c 
is the Coulomb loss rate and Fh is the hadronic loss rate 
(Paper III). In Paper III, we estimated that 



Vcm A ) 
F h = Xh ( " e „ ) 



erg cm 



erg cm - 



-1 -3 

erg s cm 



erg s cm 



(9) 



(10) 



where Xc = -7.3 x 10~ 16 and Xh = -1.5 x 10~ 17 . In these 
estimations, we assumed that the CR momentum spectrum 
is given by a power low (oc p x ), and the index is x = 3. We 
chose x — 3 because it is consistent with observations of ra- 
dio mini-halos in clusters (Paper III). Moreover, we assumed 
that the minimum momentum is p m i n c = 137 MeV, at which 
the effect of Coulomb collision is maximum (Paper III). For 
these values, we have confirmed that U C o\\ <md ri oss can be 
ignored for the results in the following sections. We note that 
the index could be as large as x = 3.5 to be consistent with 
the observations when v a t = va (Paper III). We estimated 
r c /(n e e c ) and Fh/{n e e c ) for x — 3.5 and p m i n c = 137 MeV, 
and found that X c = -7.6 x 10~ 16 and Xh = -4.9 x 10 -18 . 
Thus, H co n/(n e e c ) and F\ OBB /(n e e c ) are not much different 
from those when x — 3 and = 137 MeV. 

We also changed p m i n for x — 3. When p m i n c = 43 MeV, 



-2.7xl0~ ib and X h = -1.6xlO _iB , and whenp„ 



440 MeV, xc = -3.8 x 10 -16 and Xh = -1.2 x 10' 



mean that H co \\/(n e e c ) and ri ss/(n e e c ) are smaller than 
those when p m i n c = 137 MeV. Thus, we can ignore H co \\ 
and Fioss- It is to be noted that when x > 3, one can obtain 
r c /(n c e c ) oc Pmin regardless of x in the limit of small p m in 
(see equations [7] and [8] in Paper III). 

2.4 Cluster 

For the gravitational mass profi l e, we adopt the NFW 
model ijNavarro. Frenk. fc Whitel Il997f ). although there 
is a debate about the s lope of the central cusp (e.g. 
iFukushige fc Makinolll997l ). For the NFW profile, the mass 
distribution is written as 



M(r) = M 200 



ln(l + r/r a ) — r/r s /(l + r/r 3 ) 
ln(l + C200) - C200/ (1 + C200) 



(11) 



For this equation, we define r q as the cluster radius within 
which the average mass density is q times the critical density 
of the universe p C i(z) at redshift z. Moreover, M q is the mass 
of the cluster within r q , r s is the characteristic radius, and 
c i = r il T a- In equation (|11[1 . we set q — 200. From the 
definition, we have 

1/3 



3M„ 



4-Kqpcr(z) 



(12) 



The mass M200 is derived from the ICM tempera- 
ture outsid e the core, T out . Based on a statistical study, 
IChen et all (|2007T l obtained a scaling relation of 



M 500 = 

If we 
verted 



2.6 x 10 i4 /T 



1 / T on t \ 
UkeV/ 



Mr. 



AkeV. 

assume the NFW profile, 
to M200 using a relation of M oc 



(13) 



M: 



500 can 



be con- 

^-0.266 

l|Horner. Mushotzkv. fc Scharj|T999r ). Theoretically, C200 is 
expected to be a weakly decreasing function of M200 (e.g. 
Duffy et al.1 [20081) . However , observations hav e shown that 
it ha s a very large scatter l|Okabe et al.1 12010| ; lEttori et al.1 
2010), which may reflect a broad rang e of the formation 
epoch s of clusters with a give mass (e.g. iFuiita fc Takaharal 
1999). Thus, we fix it at C200 = 5. Using equations (|11[| - 
(|13|l . we can determine M(r) of a cluster with given T out 
and z. The gravitational acceleration for the NFW profile is 
ffNFw(r) = ~GM(r)/r 2 . 

In addition to <?nfw, we include the gravitation from 
the central cD galaxy. We use t he o ne obtained by 
iMathews. F altenba cher. fc Brighentil (|2006T ) and used in Pa- 
pers I-III: 

an -l/s 



3.206 x 10" 7 



+ 



1.861 x 10" 6 



(14) 



These 



in cgs with s = 0.9 and r in kpc. We assume that g c o does 
not depend on host clusters. Thus, the total gravitational 
acceleration in a cluster is g — jnfw + JcD . The inclusion of 
g C D does not qualitatively change the results. 



3 RESULTS 

3.1 Steady State Solutions 

First, we derive steady state solutions (d/dt — 0). The 
pro cedure is basically the same as that in previous stud - 
ies (|Kim fc Naravanl 120031 ; iGuo. Oh. fc Ruszkowskil l2008l l. 
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From equation {T}, the mass flow rate is given by M = 
Aixr 2 pu. Other three equations ©-Q can be rewritten as 



GM{r) 



r 2 dr 
r 2 dr 



- n e K(T) + H 8t , 



S c 



(15) 



(16) 



(17) 



1 d . 2 2 
P U ) 



1 d . 2 v 

^* (r ue9) = 

1 . 2- v 

We solve these ordinary differential equations for n n < r < 
r out , where n n = 3 kpc and r out = 1 Mpc using Mathemat- 
ica They can be solved as a boundary value problem. 
We impose the following boundary conditions: 

n e [r- m ) = n , (18) 

T(r in ) = T in , (19) 

T(r out ) = r out , (20) 

Pc(r in ) = P c0 . (21) 

Equations ()15p - (|17p form an eigenvalue problem in which 
M is the eigenvalue. 

When v s t = va = B/^/iirp, we need to specify 
magnetic fields B in the ICM. We assume that B = 
Bo (n e /0.016 cm -3 )'- We could not often find steady state 
solutions for too small Bo and/or too large b. Thus, we as- 
sume that B — 1.0 x 10~ 5 pG and b = 0.4, although we 
assumed that Bo = 1.0 x 10~ 5 pG and b = 2/3 in Paper III. 
The fairly small value of b may mean that the coupling be- 
tween magnetic fields and the ICM is weak, which may be re- 
alised when the magnetic fields are rather radially extended 
in a cluster. The small b is also favourable to suppress the 
development of local instabilities (see § I3.2.2p . 

We construct steady state solutions for four clusters 
with various temperatures (A17 95, A2199, A2052, and 
A2597 ), which were studied by iGuo. Oh. fc Ruszkowskil 
(2008). In Table [I] we give the boundary values no, Tin, 
and Tout when v a t = va, which were chosen to be almost 
consistent with observations (Figs. . Although there are 
no direct observations of P c o, the values of Pco/PgO, where 
P g o — Pg(ri n ), are chosen so that dP c o/dr is not positive 
and close to zero at r = r^. 

Dashed lines in Figs. HHH show the steady state solu- 
tions we obtained. Derived mass flow rates, M , are shown 
in Table [T] It is to be noted that M is the mass flow that 
passes the inner boundary (ri n = 3 kpc), and that not all 
the mass needs to fall into the central black hole. Most 
of the gas may become cold gas or be consumed by star 
formation in the central g alaxy (e.g. iBregman et al.ll200r3 ; 
iMcNamara fc Nulsenl2012h . For comparison, we show obser- 
vations of the four clusters in the figures. The steady state 
solutions can generally reproduce the observations. Since our 
model is rather simple (we do not include additional mechan- 
ical heating, for example), we think that it would be useless 
to perfectly fit the solutions with the observations. 

Fig. [5] shows the profiles of P c , Pb, u, H Bt , and S c for 
the steady state solution of A1795. The results are qualita- 
tively the same for the other clusters. The pressures (P c and 
Pb) and the infall velocity (— u) increase toward the cluster 

1 http://www.wolfram.com/ 




1(T 

r (kpc) 



Figure 1. (a) Temperature and (b) density profiles for A1795. 
Dashed lines show the steady state solution or the initial profiles 
for the numerical simulation. Solid and Dotted lines are the re- 
sults of numerical simulation at t = 40 Gyr a nd 100 Gyr, respec - 
tively. Filled circles represent observations <Ettori et alj|2002ft . 
Error bars are omitted. 



centre (Fig. [5^). Both P c and Pb are smaller than P g in 
the whole region. Although Pb > P c for r > 200 kpc, Pb 
is dynamically unimportant there. Fig. [5Jd shows that the 
slope of H st is more gentle than that of S c at r > 4 kpc. 
This reflects that some of the CRs injected in the inner core 
region stream in the ICM and heat the ICM in the outer 
core region. 



3.2 Stability Analysis 

3.2.1 Lagrangtan Perturbation Analysis 

We investigate the stability of the steady state solutions 
obtained in § 13.11 In this subsection, we study the stabil- 
ity by a Lagrangian perturbation analysis. We focus on 
the global thermal instability in the ICM. The ana lysing 
method is based on that in iKim fc Naravanl (|2003T ) and 
IGuo. Oh. fc Ruszkowskil <|2008h . 

The relation between a Lagrangian perturbation oper- 
ator A and an Eulerian perturbation 5 is 



A = 6 + £- V, 



(22) 



whe re $ is the Lagrangian disp lacement of a fluid element 
fsee IShapiro fc Teukolskvl Il983l . p. 127-147). We consider 
only radial perturbations and we define the radial compo- 
nent of £ as £ = Ar. In this case, the Lagrangian perturba- 
tion has commutation relations of 



d 
d 



d 

W 
d 



dr dr 



A 



d£d_ 

dr dr 



(23) 



(24) 
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Table 1. Cluster Parameters. 



Cluster 


z 


Vat 


Tin 


Tout 


no 


Pco/PgO 


M 








(keV) 


(keV) 


(cm -3 ) 




(Moyr- 1 ) 


A1795 


0.0632 


VA 


1.1 


7.5 


0.25 


0.15 


-27.1 






Cs 


1.5 


7.5 


0.17 


0.03 


-19.7 


A2199 


0.0309 


VA 


1.3 


5 


0.1 


0.10 


-9.2 






Cs 


1.5 


5 


0.1 


0.03 


-10.4 


A2052 


0.03549 


VA 


1.2 


4.5 


0.1 


0.10 


-8.0 






Cs 


1.3 


4.5 


0.1 


0.03 


-7.4 


A2597 


0.083 


VA 


1.3 


4.5 


0.1 


0.20 


-10.4 






Cs 


1.4 


4.5 


0.1 


0.03 


-9.5 




10' 

r (kpc) 




10' 

r (kpc) 



Figure 2. S ame as Fig. [T]but for A 2199. Filled circles represent 
observations I Johnstone et al"]|2002h . 



Figure 3. S ame as Fig. [T] but for A2052. Filled circles represent 
observations llBlanton et al.lfeplll) . 



i|Shapiro fc Teukolskv|[l983l) . 

Equations (J2J— (J4j) with k = D — H co n = ri oss 
8Pb /dr — can be rewritten as 



du 

p dt = 



dr 



+ P9 ; 

1 dPg _ 7g Pg dp 
7 9 — 1 dt 7 9 — 1 p dt 

1 dP n 



H Bt - pC . 



(25) 
(26) 



7c — 1 dt 



7c 



7 C P c dp 7c p x a fQ7 s 

H 7"c(V • WstJ = i>c , (27) 



1 p dt 7 C 



where d/cte is the Lagrangian time derivative, and v s t is 
the streaming velocity including the direction. Equation {TJ 
gives the mass flow rate, M = 4irr 2 pu. 

We linearize equations (|25|) - (|27l) . From these equations, 
and useful relations with gas density and pressure 



(28) 
(29) 



Ap = -pV ■ £ , 
AT 

APg=Pg— ~P 9 V-£ 

|Shapiro fc Teukolskvl[l983l ; iGuo. Oh. fc Ruszkowskill2008l ) , 
© 0000 RAS, MNRAS 000, 000-000 



we obtain following equations: 

^i = ^( ) 

dt 2 p dr y ^ } 



ld_ 

p dr 



7 9 P g dp\AT 



■Jg-ldt 



+ pTCr 



1 <ZP„ 



7 9 — 1 dt 7 9 — 1 p dt 



T 



+ P 



d_ 

dt 



P Z C P - Hst (V ■ - AH st = . 



d (AP C \ AS C S C AP C 

where £t = dC/dT\ p and £ p = dC/dp\r- 

Since w st = ua = B/^-inp oc p a , where cv = 6 
perturbations of u st can be represented by 

A««t = OLV st , 



(31) 
(32) 

- 0.5, 

(33) 
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1(T 

r (kpc) 



Figure 4. S ame as Fig. [T]but for A 2597. Filled circles represent 
observations ijMcNamara et al1l200lh . 



respectively. In the above equations, the density perturba- 
tion can be replaced by Ap/p — — V ■ £ (eq. [28]). 

Equations J6| and show that the source of CRs can 
be rewritten as S c = Q(r)M. We can write the perturbation 
of the source as 



AS C 



AQ AM- m 
~~Q~ M 



S c 



(36) 



where AMj„ is the perturbation of the mass flow rate at 
r — r- m . The perturbations in equation (|36|1 are 

dQ, 



AQ 



AMm 



dr 



M_ d£ 

Ui n dt 



(37) 
(38) 



where u- ln — u(ri n ). 

We take £, AT, and AP C as independent variables. We 



seek linear eigenmodes that behave as ~ e' 7 with time. 
Equations (|30l l (I32|) may be rewritten as 

, 4 / . ,/,; 1 



dr 



+ 2au + u 



dr p dr 



dt 
dr 



(39) 




10 

r(kpc) 

Figure 5. Profiles of (a) CR pressure (solid), magnetic pressure 
(dashed), and gas velocity (dot-dashed), (b) heating rate of CR 
streaming (solid) and CR injection rate (dashed) for A1795. 
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, Ap aA dP c 
AH st = -v st | a— t. - — - 
p or I or 



(34) 



-Vst 



d_ f AP 



' dr 



+ ■ 



APcdPc 
P c dr 



(35) 



-^-+pTC T + 



dP a 



7s - 1 dr 



7s 



Ja^P^dp \ AT 
1 p dr J T 



dr y s; 7s 



H st )(V-C)- Affrt 
P u d / AT\ 



1 dr V T / 



= 0, 



(40) 



_ AP e 



t— + ( <'., + u) J; (^-) + a 7c (V ■ £) + -,^.;V ■ I) 



(41) 



+7cA(V ■ u 8t ) = — 



We omit the dependence of e a hereafter. Equations 
(|41|l are first-order differential equations for the four vari- 
ables £/r, AT/T, AP C /P C , and d(£/r)/dr. We numeri- 
cally solve these equations as an eige nvalue problem, where 
the eigenvalue is the growth r ate a (|Kim fc Naravar]|2003l ; 
iGuo. Oh. fc RuszkowskH 120081 ). using Mathematical 9. Fol- 
lowing the previous studies, we set five boundary conditions. 
At the inner boundary (r = r- ln ), we give three conditions: 

e/r = i 
d 



dr \rj 



0. 



(42) 
(43) 



A(r 2 «P c ) = 2(u + « st )P c rC 

+ (a£ + u-^ + Avst) Per 2 + r 2 (u + v st )AP c 

= . (44) 

Equation ()42p is a normalisation condition, and equa- 
tion (|43l) guarantees the regularity of the solutions. Equa- 
tion (|44[) means that the perturbed CR flux is zero at the 
cluster centre. At the outer boundary (r — r out ), we adopt 
the two conditions: 



(45) 
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Figure 6. Eigcnfunctions of the stable mode for A1795 (dotted), 
A2199 (short-dashed), A2057 (long-dashed), and A2597 (dot- 
dashed) plotted as functions of radius, (a) £/r, (b) AT/T. The 
streaming velocity is v B t = va ■ 



AT = 0. 



(46) 



because the cooling time of the ICM is much larger than the 
cluster age. 

The background ICM and CR profiles are given by the 
steady state solutions derived in § I3.ll We search a that 
satisfies equations (|39 j) — (141 1) and boundary conditions l)42ft 
([46]) in the range of (10 4 Gyr)" 1 < a < (1CT 4 Gyr)" 1 . 
Since the equations are rather complicated, it takes a long 
time to find solutions. Thus, we limited ou r search to real a 
and d id n ot search imaginary a. Fo l lowin g iKim fc Naravanl 
(|2003h and lGuo. Oh. fc Ruszkowskil (|200ci ) , we first fix a and 
give arbitrarily AT/T at r = r ln . We then integrate equa- 
tions (|39|) (|41[) from r — ri n to r — r out and check whether 
the first outer boundary condition (|45fl are satisfied. If not, 
we try another AT. We use the second outer boundary con- 
dition (|46|) as a discriminant for solutions. 

As a result, we could not find any solutions with positive 
a in the above range for the four clusters. This means that 
perturbations do not grow and the steady state solutions 
are quite stable even if we do not include thermal conduc- 
tion, which has often been used to stabilise heating. On the 
other hand, we found decaying modes (a < 0). Fig. [6] shows 
the eigenfunctions for the lowest order mode. For the mode 
shown in the figure, —a ~ 0.3 Gyr -1 for the four clusters. 



3.2.2 Numerical Simulations 



In the analysis in § 13.2.11 we did not perform a complete pa- 
rameter search; we did not examine imaginary a and local in- 
stabilities, for example. Thus, in this subsection, we supple- 
mentarily study the stability of the steady state solutions us- 
ing numerical simulations. We solve equations (HJ-© with 
k = D — H co \\ = ri oss = dPs/dr = 0. The hydrodynarnic 



part of the equations is solved by a second- order advection 
upstre am splitting method (AUSM) based on lLiou fc Steffenl 
(1993). We use 300 unequally spaced meshes in the radial 
coordinate to cover a region with a radius of r out = 1 Mpc. 
The inner boundary is set at r- m = 3 kpc. The innermost 
mesh has a width of 90 pc, and the width of the outermost 
mesh is 17 kpc. While variables have zero gradients at the 
inner boundary, density and pressure are fixed at their ini- 
tial values at the outer boundary. We use the steady state 
solutions obtained in § 13. li as initial conditions (t = 0). 

Figs. lH4l show the results of the calculations. For all the 
four clusters, the ICM is stably heated at least within the age 
of the Universe (t ~ 14 Gyr). For A2199, A2057, and A2597, 
profiles are almost identical to the initial ones at t < 40 Gyr. 
However, local instabilities start developing at t > 40 Gyr 
at r ~ 50 kpc and it later affects the inner profiles. In fact, 
it has been indica ted that the heating via CR streaming is 
not locally stabl e (|Loewenstein. Zweibel. fc Begelmanlll99ll ; 
|Pfrommerll2013l ). For A1795, the local instabilities do not 
develop until t = 100 Gyr. 

Although the local instabilities develop for the three 
clusters, they are not very radical. We found that our as- 
sumption of smaller b contributes to the suppression of rapid 
development of the local instabilities. Since we assumed that 
b — 0.4, the Alfven velocity is va = B/yJAnp oc p -0 ' 1 . If 
excessive cooling increases p at the cluster centre for ex- 
ample, it reduces va there. The smaller va prevents CRs 
from escaping from the cluster centre and increases P c and 
\dP c /dr\ around the centre. Since // s t oc \dP c /dr\, the clus- 
ter centre is well-heated. The same mechanism should work 
when Vat = c s , because excessive cooling reduces the stream- 
ing velocity. Using numerical simulations, we found that the 
local instabilities we discussed here are suppressed by mod- 
erate thermal conduction of the level of ~ 1% of the Spitzer 
conductivity. 



4 DISCUSSION 

4.1 Evolution with time 

Because of possible cluster mergers and the change of AGN 
activities, it is likely that the ICM profiles sometimes deviate 
from the steady state solutions we studied above. Thus, we 
study the time scale in which the perturbed ICM profiles 
return to the steady ones. 

As an example, we choose the ICM profile of the steady 
state solution for A1795 when v st — v a as a fiducial pro- 
file (Table [T] and Fig. [1]). In this model, the efficiency e in 
equation Q was e = 2.5 x 10~ 4 . We change the efficiency 
and construct new steady state solutions, while fixing the 
boundary parameters (no, T n , Tout, and P c o). We obtain 
M = 19.0 M© yr" 1 and 39.4 M© yr" 1 for e = 5 x 10~ 4 and 
1 x 10~ 4 , respectively. Their temperature and density pro- 
files are not much different from those in Fig. [T] because the 
boundary conditions are the same. Using those profiles as 
the initial ones, we numerically solve equations HJ-((4} with 
k = D = // co ii = ri oss = dP B /dr = for v s t = va and 
e = 2.5 x 10~ 4 . This may be the situation where the heating 
efficiency of the central AGN suddenly changes at t — 0. 

In Fig. [7] we show the evolution of M, which is pro- 
portional to -//agn - The mass flow rates oscillate on a time 
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Figure 7. Evolution of M for e = 5 X 10 -4 (solid) and 1 X 10 4 
(dashed) at t = 0. 



scale of ~ 3 Gyr, which is comparable to the time scale 
of the decaying mode (— l/cr) that we studied in § 13.2.11 
The amplitude of the oscillation gradually decreases and M 
converges to that of the fiducial model of e = 2.5 x 10~ 4 
{M — —27.1 Mq yr~ ). Since the time scale of the first 
large oscillation (~ 3 Gyr) is smaller than the cluster life- 
time (~10 Gyr), we think that the discussions based on the 
steady state solutions are justified. 



4.2 Radio mini-halos 

In cool cores in some clusters, diffuse radio emissions called 
mini-halos have been observed. Our model predicts that 
the diffuse emissions originate from the CR protons heat- 
ing the cluster cores (Paper III). In particular, the observed 
radio profiles can be nicely reproduced by our model. Re- 
cently, it has also been proposed that the CR electrons 
responsible for the radio emission have been accelerated 
by turb ulence generated by the sloshing of the cool cor e 
gas fe.g. lMazzotta fe Giacintuccill2008l : IZuHone et al.ll2013h . 
Contrary to these models, our model does not require the 
sloshing because the CR protons are accelerated by the cen- 
tral AGNs. 

While there is no direct relationship between the 
mini-halos and the sloshing in our model, an ap- 
parent relationship could be observed. In our model, 
mini-halos are associated with the AGNs that are ac- 
tive in cool cores. Meanwhile, structures related to 
the sloshing (e.g. c old fronts) are often produced 
around the cool cores ifFuiita. Matsumoto. fc Wadal |2004| ; 
lAscasibar fc Markevitcri 



2006). Moreover, in our model, the 



CR electrons responsible for the radio emission are created 
through the interaction between the CR protons and the tar- 
get gas protons. Thus, more electrons are created in denser 
gas, and brighter radio emission would be observed in denser 
cores or on the denser side of a cold front associated with 
the sloshing. Of course, our model predicts that even clus- 
ters without sloshing can have the mini-halos. In the sloshing 
model, the morphology of the simulated mini-halos is com- 
plicated and depends on observing frequencies, because the 
sloshing is temporary and the CR acce leration occurs in th e 
region where turbulence is developing (|ZuHone et al.ll2013l ). 
Thus, the radio emission could be strong even not at the 
cluster centre. This is not generally true in our model, be- 



cause the radio emission tends to be strong at the cluster 
centre where the gas and CR densities are the highest. 

On the other hand, we do not think that our model 
can explain Mpc-scale radio halos that are often ob- 
served in merging clusters without cool cores, because 
AGNs are active in cool cores, and dense gas in cool 
cores is required to produce the CR electrons. Those 
Mpc-scale radio halos may be produced by electrons 
accel e rated in t urbul e nce in the ICM dBrunetti et all 
| 200ll; IPetrosianl I200l1; lOhno Takizawa, fc Shibatal |2002| 



20031 ; iBrunetti fc Lazarianl 



Fujita. Takizawa. fc Sarazin 
boilj). 



5 CONCLUSIONS 

We have investigated the stability of CR heating in clusters. 
The CRs are accelerated at the central AGNs, and heat the 
ICM through CR streaming. First, we obtained steady state 
solutions of the ICM and CR profiles so that they are con- 
sistent with observations. The solutions are obtained when 
the magnetic fields are strong enough and their association 
with the ICM is relatively weak. Then, we analysed the sta- 
bility of the solutions analytically and numerically. For the 
analytic approach, we adopted a Lagrangian stability anal- 
ysis and found that there are no globally unstable modes. 
For the numerical approach, we followed the evolution of 
the solutions for 100 Gyr and confirmed that the solutions 
are quite stable. These results, as well as the consistency 
with radio observations (Paper III), make the CR heating 
an attractive solution of the cooling flow problem. 
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APPENDIX A: SOUND VELOCITY CASE 

Here, we consider the case where the streaming velocity is 
the sound velocity of the ICM (« st = c s ). We set Pb = 
because v s t does not depend on magnetic fields. 



Al Steady State Solutions 

The steady state solutions can be obtained by replacing va 
by c 3 in § 13.11 We solve equations H15|) - (|17p with boundary 
conditions (|18|) - (|21l) . In order to match the solutions with 
observations, we adjust no, Tin, and P c o, while T out is un- 
charged from that in § 13.11 We found that P c o/P 9 o ~ 0.01- 
0.1 from comparison with radio observations (Fig. 14 in Pa- 
per III). Thus, we assume that P c o/PgO = 0.03 for all the 
four clusters. 

The dashed lines in Figs. IAHA4I show the steady state 
solutions for the four clusters. The temperature and density 
profiles generally reproduce the observations. The parame- 
ters we adopted are shown in Table [T] (v st = c s ). Fig. IA5I 
shows the profiles of P c , u, H st , and S c for the steady state 
solution of A1795. The results are qualitatively the same 
for the other clusters. Compared to Fig. [5^, P c is smaller 
because v s t = c s is larger than va, and it enhances the es- 
cape of the CRs fFig. IA5W ). Fig. lA5b shows that the heating 
(-ffst) is more widely distributed than the CR injection (S c ) 
as is the case of v s t = va (Fig. [5p). 



A2 Stability Analysis 

A2.1 Lagrangian Perturbation Analysis 



As we did in § 13.2.11 we solve equations (|39p - (|41fl under 
the five boundary conditions (|42|) (146[) . However, the terms 
including v st (equations. [33]-[3S]) must be modified. 

Since we assumed that v s t = c s oc T 1 ^ 2 , perturbations 
of v s t can be represented by 



A v at AT 



(Al) 



. fr—j \ /2£ 1 d£dT 



+V s t- 



AT n J_dT\ v^d_ / AT\ 



T V r 1 4T dr J 2 dr V T J 
The perturbation of the heating by the CR streaming is 
d£\ dP t 



(A2) 



,'1AT 

AH st = -v st I - — 



dr 



'Vst 



■ d /ap c 

c dr { P c 



+ ■ 



dr 
AP C dP c 



(A3) 



P c dr 

The background ICM and CR profiles are given by the 
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Figure Al. (a) Temperature and (b) density profiles for A1795. 
Dashed lines show the steady state solution or the initial pro- 
files for the numerical simulation. Dotted lines are the results of 
nume rical simulation at 100 Gyr. Filled circles represent observa- 
tions iEttori et alj|20o3 ). Error bars are omitted. 



steady state solutions derived in § I All We search a that 
satisfies equations (I39|) - (|41[) and boundary conditions (1-421) 
(0 in the range of (10 4 Gyr)" 1 < a < (10 -4 Gyr)" 1 . As is 
the case of v s t ~ va, we could not find any solutions with 
positive a in the above range for the four clusters. On the 
other hand, we found decaying modes [a < 0). Fig. lA6l shows 
the eigenfunctions for the lowest order mode. For the mode 
shown in the figure, —a = 0.4, 0.4, 0.2, and 0.3 Gyr -1 for 
A1795, A2199, A2057, and A2597, respectively. 

A2.2 Numerical Simulations 



As we did in § 13.2.21 we numerically solve equations ((TJl-Q 
with k = D = //coii = ri oss = dPsjdr = for v st = c s . The 
initial (t — 0) profiles are given by the steady state solutions 
we obtained in § I Al I 

Dotted lines in Figs. IA1HA4I show the profiles at t = 
100 Gyr. They are almost identical to those at t — (dashed 
lines). Contrary to the case of v s t = va (§ I3.2.2[l . local in- 
stabilities do not develop for t < 100 Gyr for all the four 
clusters. The results indicate that the steady solutions are 
very stable. 
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Figure A2. Same as Fig. lAll but for A2 199. Filled circles repre- 
sent observations l| Johnstone et al"1l2002|y 
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Figure A3. Same as Fig. lAll but for A2052. Filled circles repre- 
sent observations llBranton et alj|201lf) . 
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Figure A4. Same as Fig. lAll b ut for A25 97. Filled circles repre- 
sent observations (McNamara et al. 200ll ). 
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Figure A5. Profiles of (a) CR pressure (solid) and gas velocity 
(dot-dashed), (b) heating rate of CR streaming (solid) and CR 
injection rate (dashed) for A1795. 
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